clear
use "${data}IndirectSurveyDataReplication.dta", clear

************** COLLAPSE INDICES TO COMMUNITY LEVEL AND CHECK ACCURACY ****************

collapse (mean) PIB* IPW* PCA* poor*, by(ejnclnmid)
keep if poor_IPW_HH ~=.
******************* SIMPLE SUMMARY TABLE OF INDICES *****************************
preserve
reg PIB_leader PIB_HH
keep if e(sample)

foreach x in PIB IPW PCA {
rename `x'_leader `x'_temp 
}
estpost sum PIB_temp IPW_temp PCA_temp
estimates store A

foreach x in PIB IPW PCA {
rename `x'_temp `x'_leader
rename `x'_HH `x'_temp
}

estpost sum PIB_temp IPW_temp PCA_temp
estimates store B

foreach x in PIB IPW PCA {
rename `x'_temp `x'_HH
}

* Normalized difference (Imbens and Wooldridge 2007)

foreach x in  PIB IPW PCA {
ttest `x'_leader = `x'_HH
gen `x'_temp = (r(mu_2) - r(mu_1))/(sqrt(r(sd_2)^2 + r(sd_1)^2))
}

estpost sum  PIB_temp IPW_temp PCA_temp
estimates store C
foreach x in PIB IPW PCA { 
corr `x'_leader `x'_HH
replace `x'_temp = r(rho)
}

estpost sum  PIB_temp IPW_temp PCA_temp
estimates store D
foreach x in  PIB IPW PCA {

egen `x'sd = sd(`x'_leader)
egen `x'mu = mean(`x'_leader)
gen z`x'_leader = (`x'_leader-`x'mu)/`x'sd
drop *sd *mu
egen `x'sd = sd(`x'_HH)
egen `x'mu = mean(`x'_HH)
gen z`x'_HH = (`x'_HH-`x'mu)/`x'sd
	
reg z`x'_HH z`x'_leader, cluster(ejnclnmid) 
reg `x'_HH `x'_leader, cluster(ejnclnmid) 
replace `x'_temp = _b[`x'_leader]

}	

estpost sum PIB_temp IPW_temp PCA_temp 
estimates store E

la var PIB_temp "Simple index"
la var IPW_temp "Inverse proportion weighted index"
la var PCA_temp "Principle components index"

esttab A B C D using "${outputs}sum_indices_avgejido.tex", replace f cells("mean(fmt(3))") label style (tex) ///
collabels(none) mtitle( "Leaders" "Households" "Norm diff" "$\rho$" )  
restore

estimates clear


******************* FIGURES *******************

twoway (lpolyci IPW_leader IPW_HH, fcolor(none) clcolor(ekblue)) (line IPW_HH IPW_HH), graphregion(color(white)) ytitle("Leader IPW index") xtitle("Household IPW index") ///
legend(order(1 2) cols(1)  label(1 "95% CIs") label(2 "Kernel regression"))  yscale(range(0 .50 2))
graph 	save "${outputs}povIPWejido.gph", replace


twoway (lpolyci PIB_leader PIB_HH, fcolor(none) clcolor(ekblue)) (line PIB_HH PIB_HH), graphregion(color(white)) ytitle("Leader binary sum index") xtitle("Household binary sum index") ///
legend(order(1 2) cols(1) label(1 "95% CIs") label(2 "Kernel regression")) 
graph 	save "${outputs}povPIBejido.gph", replace

twoway (lpolyci PCA_leader PCA_HH, fcolor(none) clcolor(ekblue)) (line PCA_HH PCA_HH), graphregion(color(white)) ytitle("Leader PCA index") xtitle("Household PCA index") ///
legend(order(1 2) cols(1) label(1 "95% CIs") label(2 "Kernel regression")) 
graph 	save "${outputs}povPCAejido.gph", replace



gr combine 	"${outputs}povPIBejido.gph" "${outputs}povIPWejido.gph" "${outputs}povPCAejido.gph", ///
			graphregion(color(white))
gr export "${outputs}poverty_indices_ejido.eps", replace

estimates clear

